Enhanced Power Factor and Ultralow Lattice Thermal Conductivity Induced High Thermoelectric Performance of BiCuTeO/BiCuSeO Superlattice

Based on the first-principles calculations, the electronic structure and transport properties of BiMChO (M=Cu and Ag, Ch=S, Se, and Te) superlattices have been studied. They are all semiconductors with indirect band gaps. The increased band gap and decreased band dispersion near the valence band maximum (VBM) lead to the lowest electrical conductivity and the lowest power factor for p-type BiAgSeO/BiCuSeO. The band gap value of BiCuTeO/BiCuSeO decreases because of the up-shifted Fermi level of BiCuTeO compared with BiCuSeO, which would lead to relatively high electrical conductivity. The converged bands near VBM can produce a large effective mass of density of states (DOS) without explicitly reducing the mobility µ for p-type BiCuTeO/BiCuSeO, which means a relatively large Seebeck coefficient. Therefore, the power factor increases by 15% compared with BiCuSeO. The up-shifted Fermi level leading to the band structure near VBM is dominated by BiCuTeO for the BiCuTeO/BiCuSeO superlattice. The similar crystal structures bring out the converged bands near VBM along the high symmetry points Γ-X and Z-R. Further studies show that BiCuTeO/BiCuSeO possesses the lowest lattice thermal conductivity among all the superlattices. These result in the ZT value of p-type BiCuTeO/BiCuSeO increasing by over 2 times compared with BiCuSeO at 700 K.


Introduction
Thermoelectric (TE) materials, which can realize the direct conversion between heat and electricity, are of great significance for sustainable development. The dimensionless figure of merit ZT = S 2 σT/k is used to evaluate the performance of TE materials, where S, σ, T, and k are the Seebeck coefficient, electrical conductivity, absolute temperature, and thermal conductivity, respectively [1]. Thermal conductivity mainly consists of lattice k l and electronic thermal conductivity k e . The coupling relationship that exists between S and σ makes it difficult to tune them independently [2]. The increased electrical conductivity σ would lead to a decrease in the Seebeck coefficient [3]. Therefore, the power factor (S 2 σ) is used to evaluate the electrical transport properties [4]. According to the Wiedemann-Franc law, electrical thermal conductivity and electrical conductivity are also closely related [5,6].
Some oxide ceramic samples have been found to show high electronic transport properties [7,8], which have attracted a lot of attention. Therefore, the layered quaternary oxychalcogenide BiCuSeO has been widely investigated due to its low intrinsic thermal conductivity and excellent electrical transport properties [9,10]. The high phonon-phonon scattering efficiency of BiCuSeO leads to its extremely low lattice thermal conductivity [11,12]. In the past decades, many studies have been performed to enhance the TE performance of BiCuSeO. Li et al. introduced vacancies at Cu and Bi sites. They found that when the

Computational Methods
The superlattice structures of BiMChO (M = Ag and Cu, Ch=S, Se, and Te) were built along the c-axis. The Perdew-Burke-Ernzerhof (PBE) of generalized gradient approximation (GGA) was chosen as the exchanged-correlation potential to optimize the crystal structures on the Vienna ab initio simulation package (VASP) [25][26][27]. The convergence accuracy was set to 10 −8 eV. A 15 × 15 × 4 Monkhorst-Pack k-point mesh was used, and the kinetic energy cutoff was set to 450 eV The crystal structures were totally relaxed until the Hellmann-Feynman forces were less than 0.01 eV/Å. Then the WIEN 2K (version number: 21.1) package, which was based on the full-potential linearized augmented plane waves (FLAPW), was used to optimize the atomic positions, and forces as well as calculate band structures. The well-converged basis sets with R mt K max = 7.0 were used, where R mt and K max were the smallest muffin-tin (MT) radius and the maximum size of reciprocallattice vectors, respectively. Since the GGA methods usually underestimate the band gap value, the Tran-Blaha modified Becke-Johnson (TB-mBJ) semi-local exchange potential was employed, which has been confirmed to provide a more accurate band gap than standard GGA [28][29][30]. For our calculations, spin-orbit coupling (SOC) and spin polarization were not considered.
The BoltzTraP code was used to calculate the electrical transport properties, which is based on the semi-classical Boltzmann transport theory and relaxation time approximation (RTA) [31][32][33][34]. A non-shifted mesh with 5000 k-points was used for transport calculations. For the accuracy of phonon spectra, a 3 × 3 × 1 supercell was set and the number of atoms was 144 per supercell. The density functional perturbation theory (DFPT) was employed for the phonon calculations [16,19,35]. The phonon structure was obtained by the Phonopy package. At last, the ShengBTE code was used to calculate the lattice thermal conductivity [36]. The iterative method was used and the ngrid label in CONTROL file was set to 15 × 15 × 3, which could converge well by the converging test.

Crystal and Electronic Structures of BiMChO (M=Cu and Ag, Ch=S, Se and Te)
As can be seen in Figure 1, the superlattice structures were built along the c-axis. The top half of the superlattices consists of BiAgSeO, BiCuSO, and BiCuTeO, and the bottom half consists of BiCuSeO. The low rate of lattice mismatch means high physical stability. Then the lattice constants were optimized using the PBE exchange-correlation potentials. After the optimization, the lattice constants of BiCuSeO were slightly larger than the experiment values (see Table 1). This is mainly due to the lattice constants being usually overestimated by the PBE exchange-correlation potentials [37]. The larger atomic radii of Ag and Te atoms lead to the lattice constants of the superlattice structures of BiAgSeO/BiCuSeO (BASO/BCSO) and BCTO/BCSO being slightly larger than that of BiCuSeO. Figure 2 gives the band structures of the superlattices of BiMChO. The k-points of the first Brillouin zones of all the superlattices were automatically generated by the software XcrySDen (verison number: 1.6.2) along the high symmetry points of Γ-X-M-Γ-Z-R-A-Z The TB-mBJ exchange-correlation potential was employed for the calculations, which has been confirmed to provide a more accurate band gap than the PBE. The similar crystal structures result in the similar band structures in Figure 2. BiCuSeO and its superlattices are indirect band-gap semiconductors. However, the band dispersion and convergence near the Fermi level are different, which leads to different electrical transport properties. This is because the band structure near the Fermi level can greatly affect TE performance. There is no obvious band convergence occurring for p-type BiCuSeO near VBM (Figure 2a). However, band convergence has been clearly observed between high symmetry points Γ-X and between Z-R for superlattice BCTO/BCSO near VBM (Figure 2d). The converged bands near VBM would lead to a relatively high Seebeck coefficient for p-type BiCuTeO/BiCuSeO. The origin of the band convergence for BCTO/BCSO can be explained by Figure S2, which gives the schematic diagram of the crystal structure and its high symmetry points. The structures are the same along the directions of Γ-X and Z-R in Figure S2. Therefore, the band structures are similar along the Γ-X and Z-R directions. It can be further found that the band structure of BCTO/BCSO near VBM is mainly determined by BiCuTeO, as can be seen from Figures 2a,d and S1. This is because the Fermi level of BiCuTeO (7.32 eV) is higher than that of BiCuSeO (5.63 eV). However, the band structure near CBM is mainly determined by BiCuSeO, as can be seen by comparing Figures 2a,d and S1. Therefore, the increased energy level of VBM would lead to a decreased value of the band gap of BCTO/BCSO. This would increase electrical conductivity and the power factor for p-type BCTO/BCSO.   Figure 2 gives the band structures of the superlattices of BiMChO. The k-points of the first Brillouin zones of all the superlattices were automatically generated by the software XcrySDen (verison number: 1.6.2) along the high symmetry points of -------The TB-mBJ exchange-correlation potential was employed for the calculations, which has been confirmed to provide a more accurate band gap than the PBE. The similar crystal structures result in the similar band structures in Figure 2. BiCuSeO and its superlattices are indirect band-gap semiconductors. However, the band dispersion and convergence near the Fermi level are different, which leads to different electrical transport properties. This is because the band structure near the Fermi level can greatly affect TE performance. There is no obvious band convergence occurring for p-type BiCuSeO near VBM ( Figure  2a). However, band convergence has been clearly observed between high symmetry pointsand betweenfor superlattice BCTO/BCSO near VBM (Figure 2d). The converged bands near VBM would lead to a relatively high Seebeck coefficient for p-type BiCuTeO/BiCuSeO. The origin of the band convergence for BCTO/BCSO can be explained   Figure 2b demonstrates that an increased band gap and decreased band dispersion near VBM would lead to an enormous decrease in electrical conductivity for p-type BASO/BCSO compared with BiCuSeO. Therefore, BASO/BCSO possesses the worst power factor. Figure 2c shows that BiCuSO/BiCuSeO has a similar band structure as BiCuSeO near VBM. Nevertheless, the band gap value of BCSO/BCSO is higher than that of BiCuSeO. Therefore, the electrical conductivity of BCSO/BCSO is lower than that of BiCuSeO. This leads to a relatively worse power factor for BCSO/BCSO. Figure 3 shows that the band gap values of superlattices are between the two composed materials. The band gap value of BiCuTeO is the smallest compared with BiAgSeO and BiCuSO. Therefore, the band gap value of BCTO/BCSO is relatively smaller, which may mean that the electrical conductivity of BCTO/BCSO is the largest. This could result in the largest power factor for BCTO/BCSO. The largest band gap of BiAgSeO leads to the relatively large band gap of BASO/BCSO. This may give rise to the lowest electrical conductivity and the lowest power factor compared with all the other superlattices. BiCuTeO (7.32 eV) is higher than that of BiCuSeO (5.63 eV). However, the band structure near CBM is mainly determined by BiCuSeO, as can be seen by comparing Figures 2a,d and S1. Therefore, the increased energy level of VBM would lead to a decreased value of the band gap of BCTO/BCSO. This would increase electrical conductivity and the power factor for p-type BCTO/BCSO. Figure 2b demonstrates that an increased band gap and decreased band dispersion near VBM would lead to an enormous decrease in electrical conductivity for p-type BASO/BCSO compared with BiCuSeO. Therefore, BASO/BCSO possesses the worst power factor. Figure 2c shows that BiCuSO/BiCuSeO has a similar band structure as BiCuSeO near VBM. Nevertheless, the band gap value of BCSO/BCSO is higher than that of Bi-CuSeO. Therefore, the electrical conductivity of BCSO/BCSO is lower than that of Bi-CuSeO. This leads to a relatively worse power factor for BCSO/BCSO.   The electronic localization functions (ELFs) were calculated by the VASP package to analyze the bonding characteristics of superlattices of BiMChO in Figure 4. The isosurface value was set to 0.25 Å −3 . There is almost no difference in the bonding properties of Bi-CuSeO, which is located at the bottom parts of the superlattices (see Figure 4, black circles). However, the bonding characteristics of BiAgSeO, BiCuSO, and BiCuTeO are different. The covalent bonding states between Te-Bi atoms are more obvious than that of S-Bi, The electronic localization functions (ELFs) were calculated by the VASP package to analyze the bonding characteristics of superlattices of BiMChO in Figure 4. The isosurface value was set to 0.25 Å −3 . There is almost no difference in the bonding properties of BiCuSeO, which is located at the bottom parts of the superlattices (see Figure 4, black circles). However, the bonding characteristics of BiAgSeO, BiCuSO, and BiCuTeO are different. The covalent bonding states between Te-Bi atoms are more obvious than that of S-Bi, and Se-Bi atoms (see Figure 4, blue circles). This is due to the smallest electronegativity of Te atoms. The large atomic radius of the Te atom means a weak binding state with the Cu atom. The relatively weak intra-layer interaction leads to relatively strong inter-layer interaction between Te and Bi atoms. The stronger intra-layer binding state between S and Cu atoms leads to relatively weak inter-layer interaction between S and Bi atoms. The electronic localization functions (ELFs) were calculated by the VASP package to analyze the bonding characteristics of superlattices of BiMChO in Figure 4. The isosurface value was set to 0.25 Å −3 . There is almost no difference in the bonding properties of Bi-CuSeO, which is located at the bottom parts of the superlattices (see Figure 4, black circles). However, the bonding characteristics of BiAgSeO, BiCuSO, and BiCuTeO are different. The covalent bonding states between Te-Bi atoms are more obvious than that of S-Bi, and Se-Bi atoms (see Figure 4, blue circles). This is due to the smallest electronegativity of Te atoms. The large atomic radius of the Te atom means a weak binding state with the Cu atom. The relatively weak intra-layer interaction leads to relatively strong inter-layer interaction between Te and Bi atoms. The stronger intra-layer binding state between S and Cu atoms leads to relatively weak inter-layer interaction between S and Bi atoms. The high peak value of the effective density of states (DOS) near the Fermi level represents a high Seebeck coefficient [38]. Then, the DOS of superlattices of BiMChO was calculated. As can be seen in Figure 5, the dominant contribution atoms in the superlattices The high peak value of the effective density of states (DOS) near the Fermi level represents a high Seebeck coefficient [38]. Then, the DOS of superlattices of BiMChO was calculated. As can be seen in Figure 5, the dominant contribution atoms in the superlattices near the VBM and CBM are Cu and Bi atoms. The peak value of DOS near the VBM of BCSO/BCSO is very close to BiCuSeO. They have a close Seebeck coefficient. The peak value of DOS near the VBM of BASO/BCSO and BCTO/BCSO is smaller compared with BiCuSeO, which means a lower Seebeck coefficient. As for the CBM, the peak value of DOS of BCTO/BCSO is obviously smaller than the others. The Seebeck coefficient of n-type BCTO/BCSO is the smallest at the same carrier concentration.  The electrical transport properties are determined by the band structures near the Fermi level [39]. As shown in Figure 6, the band decomposed charge densities of BiM-ChO's superlattices were calculated by the VASP package to analyze the electrical conductivity. The top half of Figure 6 gives the band decomposed charge densities of BiMChO's  The electrical transport properties are determined by the band structures near the Fermi level [39]. As shown in Figure 6, the band decomposed charge densities of BiMChO's superlattices were calculated by the VASP package to analyze the electrical conductivity. The top half of Figure 6 gives the band decomposed charge densities of BiMChO's superlattices near the CBM. BiAgSeO and BiCuSO play a dominant role in the electrical conductivity in the superlattices of BASO/BCSO and BCSO/BCSO. BiCuSeO plays a dominant role in the BiCuTeO/BiCuSeO superlattice. At the same time, the most obvious conductive pathway (see Figure 6, red circles) is formed between Bi and Cu atoms in BCTO/BCSO, which indicates that the electrical conductivity of n-type BCTO/BCSO is the largest one. The bottom part of Figure 6 gives the band decomposed charge densities near the VBM. The top and bottom parts both contribute to the electrical conductivity in BCSO/BCSO and BCTO/BCSO. However, BiCuSeO plays a major role in electrical conductivity in the BASO/BCSO. This leads to relatively low electrical conductivity. All the atoms' contributions near the Fermi level are consistent with DOS in Figure 5.

Transport Properties of BiMChO (M=Cu and Ag, Ch=S, Se and Te)
The electrical transport properties of BiCuSeO and its superlattices were calculated using the BoltzTraP code, which is based on the semi-classical Boltzmann transport theory. The rigid band approximation (RBA) was used to simulate doping. As can be seen in Figure 7, the carrier concentration dependent / , Seebeck coefficient, and / under 300 K and 700 K are studied. Figure 7a-d show the carrier concentration dependent / of BiCuSeO and its superlattices under 300 and 700 K. The large band dispersion near the CBM leads to the large / for n-type BiCuSeO. Its maximum value was 6.76 × 10 19 Sm −1 s −1 under 300 K, corresponding to the carrier concentration of −1 × 10 21 cm −3 . However, for p-type / , BCTO/BCSO is the largest one in all the superlattices. This can be confirmed by the conductive pathway between Bi-Cu atoms in Figure 6, red circles. The maximum value of / was 2.98 × 10 19 Sm −1 s −1 under 300 K, corresponding to a carrier concentration of 1 × 10 21 cm −3 . Due to the coupling relationship between / and the Seebeck coefficient, the high / leads to a small Seebeck coefficient for p-type BCTO/BCSO. There is only BiCuSeO at the bottom half of the superlattice, which contributes to the electronic transport near the VBM. This results in extremely low electrical conductivity for the p-type BASO/BCSO.

Transport Properties of BiMChO (M=Cu and Ag, Ch=S, Se and Te)
The electrical transport properties of BiCuSeO and its superlattices were calculated using the BoltzTraP code, which is based on the semi-classical Boltzmann transport theory. The rigid band approximation (RBA) was used to simulate doping. As can be seen in Figure 7, the carrier concentration dependent σ/τ, Seebeck coefficient, and S 2 σ/τ under 300 K and 700 K are studied. Figure 7a-d show the carrier concentration dependent σ/τ of BiCuSeO and its superlattices under 300 and 700 K. The large band dispersion near the CBM leads to the large σ/τ for n-type BiCuSeO. Its maximum value was 6.76 × 10 19 Sm −1 s −1 under 300 K, corresponding to the carrier concentration of −1 × 10 21 cm −3 . However, for p-type σ/τ, BCTO/BCSO is the largest one in all the superlattices. This can be confirmed by the conductive pathway between Bi-Cu atoms in Figure 6, red circles. The maximum value of σ/τ was 2.98 × 10 19 Sm −1 s −1 under 300 K, corresponding to a carrier concentration of 1 × 10 21 cm −3 . Due to the coupling relationship between σ/τ and the Seebeck coefficient, the high σ/τ leads to a small Seebeck coefficient for p-type BCTO/BCSO. There is only BiCuSeO at the bottom half of the superlattice, which contributes to the electronic transport near the VBM. This results in extremely low electrical conductivity for the p-type BASO/BCSO. of DOS near the CBM for BCTO/BCSO is only 3.7 states/eV (see Figure 5). This leads to the low Seebeck coefficient for n-type BCTO/BCSO. Its Seebeck coefficient was −77 under 700 K, corresponding to a carrier concentration of −1 × 10 21 cm −3 . The high / and large Seebeck coefficient mean a large power factor. Therefore, the power factor relative to the relaxation time, / , is calculated in Figure 7i-l. Due to the smaller / and close Seebeck coefficient compared with BiCuSeO for n-type BiAgSeO/BiCuSeO and BiCuSO/BiCuSeO, their values of / are low. The / and Seebeck coefficient of n-type BCTO/BCSO are both smaller than that of n-type BiCuSeO, which results in its small value of / . As for the p-type BCTO/BCSO, its extremely higher / and slightly smaller Seebeck coefficient compared with p-type BiCuSeO lead to relatively large / . The maximum value of / of p-type BCTO/BCSO (Bi-CuSeO) was 5.7 × 10 11 WK −2 m −1 s −1 (4.9 × 10 11 WK −2 m −1 s −1 ) at 700 K, corresponding to the carrier concentration of 0.5 × 10 21 cm −3 (1 × 10 21 cm −3 ). The phonon spectra were calculated using the Phonopy and VASP packages to analyze the phonon transport properties and stability. In order to ensure the accuracy of calculations, a 3 × 3 × 1 supercell was set. As can be seen in Figure 8, there is no obvious imaginary frequency for BCSO/BCSO and BCTO/BCSO, which means high stability. However, there is more imaginary frequency around the high symmetry points for BASO/BCSO. The poor stability for BASO/BCSO limits its application. The frequency gaps for BASO/BCSO and BCTO/BCSO are 1.72 THz and 1.65 THz, respectively. Figure 7e-h show the carrier concentration-dependent Seebeck coefficient of BiCuSeO and its superlattices under 300 and 700 K. The peak value of DOS near the VBM is close to BiCuSeO and BCSO/BCSO, which was 6.9 states/eV (see Figure 5). This leads to a similar Seebeck coefficient for p-type BiCuSeO and BCSO/BCSO, both of which were 200 µVK −1 under 700 K, corresponding to a carrier concentration of 1 × 10 21 cm −3 . The low peak value of effective DOS near the VBM of BASO/BCSO and the slowly increased peak value near the VBM of BCTO/BCSO lead to their low Seebeck coefficient. The rapidly increased peak value of DOS near the Fermi level means a relatively high Seebeck coefficient. The peak value near the CBM of BiCuSeO, BASO/BCSO, and BCSO/BCSO are close, which is about 6.7 states/eV (see Figure 5). Their Seebeck coefficient was −95 µVK −1 under 700 K, corresponding to a carrier concentration of −1 × 10 21 cm −3 . The first peak value of DOS near the CBM for BCTO/BCSO is only 3.7 states/eV (see Figure 5). This leads to the low Seebeck coefficient for n-type BCTO/BCSO. Its Seebeck coefficient was −77 µVK −1 under 700 K, corresponding to a carrier concentration of −1 × 10 21 cm −3 .
The phonon spectra were calculated using the Phonopy and VASP packages to analyze the phonon transport properties and stability. In order to ensure the accuracy of calculations, a 3 × 3 × 1 supercell was set. As can be seen in Figure 8, there is no obvious imaginary frequency for BCSO/BCSO and BCTO/BCSO, which means high stability. However, there is more imaginary frequency around the high symmetry points for BASO/BCSO. The poor stability for BASO/BCSO limits its application. The frequency gaps for BASO/BCSO and BCTO/BCSO are 1.72 THz and 1.65 THz, respectively. The Grüneisen parameters and phonon group velocities of the superlattices are calculated to further understand the phonon transport properties. The large extremum of Grüneisen parameters and small group velocities means high anharmonicity [26]. The strong anharmonicity leads to the high efficiency of phonon scattering and low lattice thermal conductivity. Figure 9 shows the Grüneisen parameters dependent on the vibration frequency of phonons of the superlattices. The extremums of Grüneisen parameters of BASO/BCSO, BCSO/BCSO, and BCTO/BCSO are 16.9, 4.9, and 19.8, respectively. The large extremum of Grüneisen parameters of BCTO/BCSO means high anharmonicity. For three superlattices, the ZA branches play a dominant role in the phonon vibrations. The Grüneisen parameters and phonon group velocities of the superlattices are calculated to further understand the phonon transport properties. The large extremum of Grüneisen parameters and small group velocities means high anharmonicity [26]. The strong anharmonicity leads to the high efficiency of phonon scattering and low lattice thermal conductivity. Figure 9 shows the Grüneisen parameters dependent on the vibration frequency of phonons of the superlattices. The extremums of Grüneisen parameters of BASO/BCSO, BCSO/BCSO, and BCTO/BCSO are 16.9, 4.9, and 19.8, respectively. The large extremum of Grüneisen parameters of BCTO/BCSO means high anharmonicity. For three superlattices, the ZA branches play a dominant role in the phonon vibrations. Figure 10 shows phonon group velocities of BiMChO and its superlattices, which are dependent on the vibration frequency of phonons. The top half consists of the phonon group velocities along the a-direction. The phonon group velocities of BCTO/BCSO are 3.63 km/s and 3.68 km/s for ZA and LA branches, respectively, which is slightly smaller compared with BASO/BCSO and BCSO/BCSO. The bottom half consists of the phonon group velocities along the c-direction. The phonon group velocity of BiCuSeO is 3.43 km/s for the LA branch, which is much larger compared with that of BCTO/BCSO (1.44 km/s). The low phonon group velocities along the c-direction mean that BCTO/BCSO has low lattice thermal conductivity. strong anharmonicity leads to the high efficiency of phonon scattering and low lattice thermal conductivity. Figure 9 shows the Grüneisen parameters dependent on the vibration frequency of phonons of the superlattices. The extremums of Grüneisen parameters of BASO/BCSO, BCSO/BCSO, and BCTO/BCSO are 16.9, 4.9, and 19.8, respectively. The large extremum of Grüneisen parameters of BCTO/BCSO means high anharmonicity. For three superlattices, the ZA branches play a dominant role in the phonon vibrations.     At last, the temperature-dependent lattice thermal conductivity ( ) of BiMChO and its superlattices are calculated using the ShengBTE code, which is based on the Boltzmann transport theory. As can be seen in Figure 11, the of BiCuSeO is the largest one compared with all the superlattices. At 300 K and 800 K, the of BiCuSeO are 0.85 W/m K and 0.33 W/m K, respectively. It is very close to experiment [7]. The of the superlattices of BiMChO decreased. The of BCTO/BCSO is extremely low compared with BiCuSeO. Its is 0.21 W/mK and 0.08 W/mK at 300 K and 800 K, which is only 25% of BiCuSeO. At last, the temperature-dependent lattice thermal conductivity (k l ) of BiMChO and its superlattices are calculated using the ShengBTE code, which is based on the Boltzmann transport theory. As can be seen in Figure 11, the k l of BiCuSeO is the largest one compared with all the superlattices. At 300 K and 800 K, the k l of BiCuSeO are 0.85 W/m K and 0.33 W/m K, respectively. It is very close to experiment [7]. The k l of the superlattices of BiMChO decreased. The k l of BCTO/BCSO is extremely low compared with BiCuSeO. Its k l is 0.21 W/mK and 0.08 W/mK at 300 K and 800 K, which is only 25% of BiCuSeO.
At last, the temperature-dependent lattice thermal conductivity ( ) of BiMChO and its superlattices are calculated using the ShengBTE code, which is based on the Boltzmann transport theory. As can be seen in Figure 11, the of BiCuSeO is the largest one compared with all the superlattices. At 300 K and 800 K, the of BiCuSeO are 0.85 W/m K and 0.33 W/m K, respectively. It is very close to experiment [7]. The of the superlattices of BiMChO decreased. The of BCTO/BCSO is extremely low compared with BiCuSeO. Its is 0.21 W/mK and 0.08 W/mK at 300 K and 800 K, which is only 25% of BiCuSeO. Figure 11. The temperature-dependent lattice thermal conductivity (k l ) of BiCuSeO and its superlattices.
By the above analysis of electronic and phonon transport properties of superlattices, it can be seen that the σ/τ of p-type BCTO/BCSO increased by two times compared with BiCuSeO. Though its Seebeck coefficient decreased slightly, its S 2 σ/τ increased by 15%. At the same time, the k l of BCTO/BCSO decreased to 25% of BiCuSeO, which is only 0.09 W/mK at 700 K. We assume that the relaxation time of BCTO/BCSO is the same with BiCuSeO. The increased σ/τ of BiCuTeO/BiCuSeO leads to the increase of k e . According to the work on BiCuSeO from Wang et al. [12], the k e of BCTO/BCSO is about 0.43 W/mK at 700 K, which is twice that of BiCuSeO. The relatively low thermal conductivity and large power factor lead to its high ZT value. The maximum ZT value of p-type BCTO/BCSO was 1.3, which increased by over 2 times compared with BiCuSeO.

Conclusions
The electronic structures and transport properties of superlattices of BiMChO (M=Ag and Cu, Ch=S, Se, and Te) were studied by first-principles calculation and Boltzmann transport theory. The following aspects are summarized. (1) The decreased band gap and large band dispersion near the VBM lead to the large σ/τ for p-type BCTO/BCSO, which is two times larger compared with BiCuSeO. Though the small peak value of effective DOS near the VBM results in its low Seebeck coefficient, the S 2 σ/τ of p-type BCTO/BCSO increased by 15% compared with BiCuSeO. (2) The relatively larger Grüneisen parameter of BCTO/BCSO means high inharmonicity for phonon vibration and high efficiency for phonon scattering. At the same time, its low phonon group velocity along the c-direction also leads to low lattice thermal conductivity. The k l of BCTO/BCSO is 0.09 W/mK at 700 K, which is much smaller than that of BiCuSeO. (3) The ZT value of p-type BCTO/BCSO increases by over 2 times and arrived at 1.3 at 700 K. This work realizes the increase of power factor and decrease of lattice thermal conductivity simultaneously. The detailed analyses of electronic and phonon properties provide deep insight into the TE properties of superlattices. Though the construction of superlattice is relatively difficult in experiments compared with doping and some other common synthetic methods, we hope that our theoretical study may stimulate the research on the application of superlattice structures of BiMChO.